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ABSTRACT 

Finite  element  problems  can  often  be  partitioned  into  subproblems 
corresponding  to  subregions  into  which  the  region  is  partitioned  or  from 
which  is  originally  assembled.  In  this  paper,  iterative  methods  are  developed 
in  which  the  interactions  across  the  curves  or  surfaces  which  divide  the  region 
are  handled  by  a  conjugate  gradient  method.  Recent  progress  on  iterative 
substructuring  methods  for  elliptic  problems  in  two  dimensions  is  discussed  as 
well  as  a  connection  between  certain  of  these  methods  and  multigrid 
algorithms. 

1.  Introduction  Iterative  substructuring  methods  are  domain  decomposition  methods  in 
which  the  often  very  large  linear  systems  of  algebraic  equations,  which  arise  when  elliptic 
problems  are  discretized  by  finite  differences  or  finite  elements,  are  solved  with  the  aid  of 
solvers  for  the  same  or  similar  equations  defined  on  subregions.  The  subregions,  also  called 
substructures,  do  not  overlap.  The  interaction  between  the  substructures,  to  enforce  various 
continuity  requirements,  is  handled  by  an  iterative  method,  typically  of  conjugate  gradient 
type. 

In  this  paper,  we  will  assume  that  the  subproblems  are  solved  exactly,  but  we  note  that 
certain  inexact  solvers  are  well  understood;  cf.  Bramble,  Pasciak  and  Schatz  [6]  ,  Dryja  , 
Proskurowski  and  Widlund  [IS]  and  Section  2  of  this  paper. 

The  new  results  reported  here  have  recently  been  obtained  in  joint  work  with 
Maksymilian  Dryja,  who  in  his  report  to  this  conference  focuses  on  our  results  for  three 
dimensional  problems;  see  Dryja  [12]  .  Only  problems  in  two  dimensions  are  discussed,  in 
this  paper.  We  consider  finite  element  approximations  of  self  adjoint,  linear,  elliptic 
problems  formulated  in  the  traditional  variational  form  as  minimization  problems.  In 
developing  the  theory,  we  will  need  an  extension  theorem  for  finite  element  functions  on 

V 

regions  with  Lipschitz  boundaries,  similar  to  well  known  results  on  Sobolev  spaces;  cf.    Necas 
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[19]  or  Stein  [21].  We  know  how  to  prove  such  a  result  for  general  geometric  configurations 
and  triangulations,  but  only  for  conforming  finite  elements;  c£.  Widlund  [24].  Therefore,  for 
the  time  being,  we  confine  our  study  to  finite  element  spaces  which  are  subspaces  of  the 
Sobolev  space  appropriate  for  the  variational  form  of  the  elliptic  equation  being  considered. 

Good  iterative  substructuring  algorithms  lead  to  improved  methods  already  for  uni- 
processor systems  and  two  substructures.  In  the  context  of  parallel  computing,  we  are  of 
course  primarily  concerned  with  problems  with  many  substructures.  In  our  discussion,  two 
themes  emerges:  When  specialized  to  the  case  of  two  substructures  the  algorithms  should  be 
optimal,  i.e.  their  rates  of  convergence  should  be  bounded  independently  of  the  size  of  the 
problem.  Secondly,  we  should  also  assure  that,  in  each  step  of  the  algorithm,  there  is  enough 
global  transportation  of  information.  These  issues  are  discussed  in  Section  2,  where  we 
review  two  basic  optimal  algorithms  for  the  case  of  two  substructures;  cf.  Bj^rstad  and 
Widlund  [3,4],  Bramble,  Pasciak  and  Schatz  [6],  Dryja  [9,10]  or  Widlund  [23,24].  We  also 
consider  the  case  of  many  substructures  and  give  a  lower  bound  on  the  rate  of  convergence  of 
iterative  substructuring  methods  for  which,  in  each  iteration  step,  information  is  exchanged 
only  between  neighboring  substructures. 

In  Section  3,  we  collect  a  number  of  results  on  Sobolev  and  finite  element  spaces,  which 
are  fundamental  for  the  proofs  of  our  results  as  well  as  for  an  understanding  of  the 
algorithms  and  their  structure.  Some  of  these  tools  do  not  have  counterparts  in  the  continuous 
case.  The  central  technical  result  is  given  in  Lemma  3.2.  In  the  last  Section,  we  introduce  a 
number  of  algorithms  and  establish,  by  using  the  tools  of  Section  3,  that  they  are  almost 
optimal.  To  simplify  the  presentation,  we  develop  these  methods  and  the  theory  primarily  for 
second  order  problems,  but  as  demonstrated  by  a  discussion  of  a  biharmonic  problem,  our 
techniques  are  also  applicable  to  higher  order  equations. 

We  note  that  the  first  algorithm  and  theoretical  results  of  this  kind  are  due  to  Bramble, 
Pasciak  and  Schatz  [7].  We  present  a  different  derivation  of  their  algorithm.  Our  proofs  of 
the  almost  optimality  of  their  algorithm  and  others  are,  in  our  opinion,  shorter  and  add  new 
insight. 

For  a  further  discussion  of  the  literature  cf.  Bj^rstad  and  Widlund[4],  where  various 
algorithms  for  the  case  of  two  substructures  are  discussed  in  detail.  The  interest  in  iterative 
substructuring  algorithms  has  grown  rapidly  recently.  A  number  of  papers  have  appeared 
since  that  paper  was  given  its  final  form;  cf.  Widlund  [25],  where  further  information  on  the 
literature  is  given.  A  more  detailed  discussion  of  the  results  presented  in  this  paper  will  be 
given  in  Dryja,  Proskurowski  and  Widlund  [16]. 

We  will  not  present  any  numerical  results  in  this  paper.  It  deserves  to  be  mentioned, 
however,  that  in  our  experience  7-10  iteration  are  sufficient  to  decrease  the  error  by  a  factor 
10*  for  many  problems. 

2.  Basic  algorithms  and  a  lower  bound.  Let  H  be  a  Lipschitz  region  and  let  its 
boundary  be  the  union  of  Fq  and  Fy  on  which  essential  (Dirichlet)  and  natural  (Neumann) 
boundary  conditions  are  imposed.  A  linear,  self  adjoint,  elliptic  problem  of  order  2€  is  given 
in  variational  form  as 

an('«.v)  =/(v),  V  V  €  Vo(n). 


The  solution,  u  and,  in  the  case  of  higher  order  equations,  €-1  of  its  normal  derivatives, 
take  on  given  values  on  Fq.  The  space  VQ(ft)  is  the  subspace  of  V{Cl),  e.g.  the  Sobolev  space 
H^(il),  or  (H^(il)y,  with  zero  data  on  Fq.  The  linear  functional  /(v)  is  obtained  from  the 
right  hand  side  of  the  elliptic  equation  and  the  natural  boundary  conditions  defined  on  F  ^  by 
using  a  Green's  formula.  We  assume  that  the  bilinear  form  a^{u,v)  satisfies  the  following 
standard  conditions, 


c\\u\\l^^  aiu,u)     , 
|a„(«,v)|  S  C  ||u||u  l|v||v„. 


(2.1) 


where  c  and  C  are  positive  constants.  It  follows  immediately  that  the  discrete  problem 
obtained  by  a  Galerkin  procedure  using  a  conforming  finite  element  space  V''  C  V  ,  has  a 
positive  definite,  symmetric  stiffness  matrix  K. 

Associated  with  the  finite  dimensional  space  V''  is  a  triangulation  of  fl  which  can  be 
quite  general,  except  that  we  insist  that  it  is  regular,  i.e.  there  is  a  uniform  bound  on  h^p^. 
Here  h^  is  the  diameter  and  p;^  the  radius  of  the  largest  inscribed  circle  of  the  element  K. 
No  essential  complications  arise  if  we  allow  some  of  the  triangles  to  be  curved  as  in 
isoparametric  methods;  cf.  Bernardi  [2]  and  Ciarlet  [8].  We  could  also  develop  the  same 
theory  for  quadrilateral  elements. 

To  be  able  to  discuss  the  case  of  two  substructures,  we  partition  fl  into  nonoverlapping 
subregions  flj  and  fl,.  The  curve  which  dissects  the  region  is  Fj  =  fl,  D  flj.  the  intersection 
of  the  closures  of  fl,  and  fl,,  and  the  boundaries  of  fl,  fl,  and  flj  are  given  as 
r,  U  Fj,  F,  U  Fj  and  Fj  U  F3,  respectively.  We  always  assume  that  F3  follows  element 
boundaries  of  the  triangulation.  We  use  basis  functions  in  V''  which  satisfies  the  standard 
conditions,  e.g.  they  vanish  outside  the  triangles  with  which  their  respective  degrees  of 
freedom  are  associated;  cf.  Ciarlet  [8].  The  elements  of  K  are  obtained  as  t,y  =  a^{^,,^j), 
where  <|),  and  <|>  are  basis  functions.  This  leads  to  zero  blocks  in  the  stiffness  matrix  and  a 
linear  system  of  equations  corresponding  to  the  Galerkin  equation  which  is  of  the  form. 


Kx  = 


0 

^13' 

*1 

"1 

^22 

^23 

JCj 

= 

h 

^^23 

^33 

,^3. 

h 

(2.2) 


The  matrix  /T,,  represents  the  couplings  between  pairs  of  degrees  of  freedom  associated  with 
the  open  set  flj  as  well  as  the  intersection  of  F^y  and  fl,,  K^^  the  couplings  between  pairs  with 
one  element  in  this  set  and  the  other  associated  with  F3,  etc.  The  right  hand  sides  are 
obtained  from  /(((>;)  and  the  Dirichlet  data. 

When  considering  iterative  substructuring  algorithms,  we  assume  that  it  is  acceptable  to 
solve  discrete  problems  on  the  subregions,  with  appropriate  boundary  conditions  on  their 
boundaries.  We  can  therefore  concentrate  on  the  case  when  b^  and  fc,  are  zero,  since  we  can 
reduce  system  (2.2)  to  such  a  form  in  a  preliminary  step.  We  call  the  finite  element  functions 
corresponding  to  such  special  right  hand  sides  piecewise,  discrete  harmonic  functions.    By 


block  Gaussian  elimination,  the  system  (2.2)  then  reduces  to 

Sx^  =  (ATjj  -  ATij  K^^   K^^  -  ATjs  ^22   ^23)-*3  ^  ^3  (2-3) 

The      so      called      Schur      complement      S      is      given      by      5  =  5''^  +  5<^\      where 
J(')  =  AT^^  -  irfj  Jf  fi'  Jr,3,  is  the  the  Schur  complement  associated  with  the  problem  , 


ICiD^lU  = 


^11  ^13 


0  . 

(2.4) 


Equation  (2.4)  is  the  finite  element  approximation  of  the  elliptic  problem  restricted  to 
n,  with  a  natural  boundary  condition  added  on  Fj.  The  elements  of  K\^^  are  a^{<^^,  <^j), 
where  4>,  and  <j)  are  basis  functions  associated  with  degrees  of  freedom  on  Fj.  The  matrix 
Jfjj  in  (2.1)  is  thus  the  sum  of  /f^y  and  K'^^^ ,  the  contributions  from  ft,  and  flj  respectively.  It 
is  also  clear  that  K  can  be  derived  easily  from  AT'"  and  K'-'^K 

While  this  discussion  might  seem  somewhat  special,  since  it  refers  to  a  problem  cut  in 
two,  matrices  of  a  general  structure  as  in  (2.2)  appear  very  frequentiy  when  two  subproblems 
are  merged.  In  industrial  finite  element  work  a  multi-level  substructuring  approach  is  often 
adopted  in  which,  at  successive  levels,  pairs  of  increasingly  large  substructures  are  merged; 
cf.  Bell,  Hatlestad,  Hansteen  and  Araldsen  [1]. 

As  shown,  e.g.  in  Bj(Arstad  and  Widlund  [4]  and  Widlund  [24],  the  condition  number  of 
5  goes  to  infinity  when  the  problem  size  increases  and  a  preconditioner  for  S  should 
therefore  be  found.  Not  surprisingly,  5^"  serves  well  in  a  conjugate  gradient  algorithm  which 
solves  equation  (2.3).  As  always  the  rate  of  convergence  is  determined  by  the  stationary 
values  of  a  generalized  Rayleigh  quotient 

xlSxj/xlS^'^Xj  (2.5) 

In  particular,  the  method  is  optimal  if  the  values  of  the  Rayleigh  quotient  lie  in  a  fixed 
interval.    We  denote  the  corresponding  condition  number  by  k(5(5'")~'). 

It  is  important  to  note  that  this  Rayleigh  quotient  is  the  ratio  of  the  strain  energy  x^Kx 
which,  by  a  simple  calculation,  is  equal  to  xlSx^  for  a  piecewise,  discrete  harmonic  function 
represented  by  the  vector  x,  and  the  strain  energy  attributable  to  the  subregion  flj.  A  lower 
bound  is  thus  obtained  immediately  but  an  upper  bound  requires  an  extension  theorem;  see 
Widlund  [24]  for  a  discussion  and  a  proof  of  such  a  theorem  for  general,  conforming  finite 
element  spaces. 

Remark.  The  method  just  introduced  could  be  called  the  Neumann-Dirichlet  algorithm.  A 
Dirichlet-Neumann  preconditioner  can  also  be  introduced.  To  simplify  our  notations,  we  only 
consider  a  second  order  problem.  We  first  solve  a  problem  in  C^^  with  a  Dirichlet  condition 
on  F3.  The  function  is  then  continued  to  ftj  by  solving  a  Neumann  problem,  while  assuring 
that  the  jump  in  the  normal  derivative  across  F3  takes  on  the  right  value.  Finally  the  residual 
is  computed  as  the  difference  between  the  function  values  on  F3.  The  iteration  operator  for 
this  algorithm  turns  out  to  be  (i'''')"'5  and  there  is  thus  a  close  relationship  between  these 
algorithms. 
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For  a  discussion  of  a  number  of  non-optimal  algorithms;  see  Bj6rstad  and  Widiund  [4]. 
Among  these  are  the  rather  natural  algorithms  which  use  the  values  of  u  or  du/dn  on  Fj  as 
the  basic  unknown  in  the  iteration  and  proceed  simultaneously  to  compute  the  solution  in  the 
two  subregions  with  the  Dirchlet  or  Neumann  boundary  values. 

Before  turning  to  our  other  basic  algorithm,  we  write  down  the  first  preconditioner  in 
matrix  form.  A  vector  x,  corresponding  to  a  piecewise,  discrete  harmonic  function  is 
constructed  by  first  solving  equation  (2.4).  Then  values  of  the  degrees  of  freedom  in  flj  ^^ 
found  by  solving  the  equation 


K, 


22 


^23  -"^3- 


(2.6) 


We  can  represent  the  operations  corresponding  to  (2.3)  and  (2.5)  in  terms  of  a  linear  system 
of  equations. 


(2.7) 


^■n   0 

^.3 

^1 

■  0 

0     a:.. 

^23 

^2 

= 

0 

K\,   0 

.    ''3, 

c 

After  a  permutation,  this  matrix  resembles  a  block  Gauss-Seidel  operator,  but  we  note 
that  one  of  its  diagonal  blocks  differs  in  an  important  way  from  that  which  appear  in  that 
method.  The  splitting  we  are  considering  here  is  of  course  based  on  problems  defined  on 
substructures. 

The  matrix  in  (2.7)  is  not  symmetric,  but  as  can  be  shown  by  block  Gaussian 
elimination,  the  restriction  of  this  operator  to  the  subspace  of  vectors  which  represent 
functions  which  are  discrete  harmonic  in  H^  is  symmetric.  By  introducing  an  additional 
factor,  we  obtain  a  preconditioner  which,  by  a  simple  calculation,  is  symmetric  on  the  entire 
space, 


ii  = 


0 

/ 


0  K21K22 


0 

0 

'  / 

0     K 


0 


22    ^23 


•13 


0     AT^^J 


0     K 


0     K, 


13 
22    *^23 


^13    ^^23    ^33.* 


(2.8) 


where  Jfjj  =  K\^^  +  kI^  JPjV  ^23-  ^*  °°*^  ^^^  *^  ^'"^  factor  in  the  factored  form  of  K  , 
represents  an  operation  which  is  trivial  in  the  case  when  the  second  subvector  of  the  right 
hand  side  vanishes. 

When  this  method  is  used,  an  additional  Dirichlet  problem,  on  one  of  the  subregions, 
has  to  be  solved  in  each  iteration  but  the  advantage  is  that  exact  solvers  for  the  subproblems 
no  longer  are  required.  Thus  if  we  wish  to  solve  a  linear  system  of  equations  of  the  form 
Kx  =  b,  where 

c  x^Kx    s  x^Kx  s,     Cx^Kx, 

then  we  can  use  the  matrix  given  in  (2.8)  as  a  preconditioner.  It  is  easy  to  see  that 

k(jp"'  k)  rs  k(a"'  k)  k  (a:->  k)  ^  k(^~'a:)  c/c  . 
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The  condition  number  oi  K'^K  is  the  same  as  that  of  5"*  S.  This  can  be  seen  by  using 
the  same  block  Cholesky  factors  to  turn  K  and  K  into  block  diagonal  matrices.  The 
appropriate  block  lower  factor  is 


/ 

0   0 

/ 

0 

0 

0 

/   0 

0 

/ 

0 

i^ll 

0    / 

/ 

23     22 

/ 

The  resulting  block  diagonal  matrices  are 


^11   0       0 

0      K22   0 

and 

0      0'  s 

'Kn 

0 

0 

0 

^22 

0 

0 

0 

5(1) 

The  Rayleigh  quotient  can  then  be  written  as 

(yl^nVi  +  ^2^22^2  +  ylsy  iV(y\K,,y^  +  ylK22y2  +  y[5"V3). 

and  the  result  follows.    We  note  that  this  study  of  preconditioners  using  inexact  solvers 
requires  tools  of  linear  algebra  only. 

Similar  constructions  are  possible  for  other  preconditioners.  This  idea  was  in  fact  first 
introduced  in  Bramble,  Pasciak  and  Schatz  [6]  for  the  case  of  two  subregions  and  the  second 
standard  preconditioner  which  we  will  now  introduce.  To  simplify  our  notations,  we  will 
confine  our  discussion  to  the  case  of  a  second  order  problem  with  Dirichlet  boundary 
conditions  on  3fl  =  T,  U  Fj. 

A  discrete  harmonic  function  is  uniquely  defined  by  its  trace,  i.e.  its  boundary  values, 
and  it  is  therefore  natural  to  construct  a  preconditioner  of  S  by  using  an  appropriate  norm  of 
a  function  defined  on  F.  It  is  known  that  in  the  continuous  case  the  W"^(3n)-norm  of  the 
boundary  values  is  equivalent  to  the  //'(H)  norm  of  the  function  which  is  the  harmonic 
extension  of  these  boundary  values.  This  result  can  be  extended  to  the  discrete  case  by  using 
the  trace  theorem  and  the  extension  theorem  previously  mentioned;  cf.  Widlund  [24].  The 
//'*(3fl)-norm  can  be  defined  in  terms  of  the  square  root  of  the  differential  operator 
-{d/dsY.  In  the  case  at  hand  the  only  non-trivial  part  of  the  trace  is  confined  to  F3.  The 
corresponding  preconditioner,  denoted  by  J  in  Bj6rstad  and  Widlund  [4],  is  essentially  the 
square  root  of  the  finite  element  approximation  of  this  differential  operator  on  Fj,  using  zero 
Dirchlet  values  at  the  end  points  of  F3  and  the  restriction  of  the  finite  element  space  V*  to  Fj. 
The  related  space,  H^^iX i),  will  be  further  discussed  in  Section  3. 

It  deserves  to  be  noted  that  in  simple  cases  the  operator  J  can  be  represented  with  the 
aid  of  fast  Fourier  transforms;  cf.  Bjcirstad  and  Widlund  [4]  or  Bramble,  Pasciak  and  Schatz 
[6]. 

We  conclude  this  Section  by  considering  the  case  of  many  substructures.  Let  fl,  be  a 
substructure  of  H  and  let  AT*''  and  x^^^^  K'-'^x'-'^  be  the  stiffness  matrix  and  the  strain  energy 
corresponding  to  fl,.  We  denote  by  F,y  a  curve  that  is  the  intersection  of  fl,  and  fl^,  the 
closure  of  two  neighboring  substructures.  Without  seriously  decreasing  the  generality  of  our 
arguments,  we  will  assume  that  the  substructures  are  triangular.  We  thus  have  two  levels  of 
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triangulations  of  fl,  in  terms  of  the  substructures  and  the  finite  element  triangulation.  We 
will  refer  to  these  as  the  coarse  and  fine  triangulations.  Our  estimates  of  the  rate  of 
convergence  will  be  given  in  terms  of  log(H/h),  where  H  and  h  are  the  diameters  of  a 
substructure  and  an  element,  respectively.  Our  arguments  will  be  local;  one  substructure  and 
its  next  neighbors  will  be  examined  at  a  time.  We  can  think  of  this  logarithmic  quantity  as  a 
measure  of  the  maximum  number  of  refinement  steps  required  when  going  from  any 
substructure  to  its  finest  refinement  level  typically  by  repeatedly  cutting  triangles  in  four;  cf. 
Yserentant  [26] 

It  is  well  known  that  for  elliptic  problems,  a  residual  of  modestly  low  frequency  will 
result  in  an  error  which  cannot  be  neglected  in  any  part  of  the  region.  Therefore  any  iterative 
method  in  which  information  is  exchanged  only  between  neighboring  substructures  must 
necessarily,  for  certain  initial  errors,  require  a  number  of  steps  which  is  at  least  equal  to  the 
diameter  of  the  dual  graph.  We  recall  that  the  dual  graph  corresponding  to  a  triangulation  is 
constructed  by  introducing  a  vertex  for  each  triangle  and  an  edge  between  any  pair  of  nodes 
which  represents  triangles  which  share  an  edge  F,^  of  a  substructure.  The  diameter  of  the 
graph  is  the  maximum  distance  between  pairs  of  nodes,  where  the  distance  is  measured  by 
the  length  of  the  shortest  path  between  nodes.  In  a  case  where  the  diameter  of  the  region  is 
one,  this  diameter  is  typically  on  the  order  of  1/H.  It  is  now  easy  to  show,  by  an  argument  of 
contradiction  and  a  well  known  estimate  of  the  rate  of  convergence  of  conjugate  gradient 
algorithm,  that  the  condition  number  of  the  iteration  operator  is  at  least  on  the  order  of  l/H^. 
It  is  shown,  in  Section  4,  that  it  is  easy  to  obtain  algorithms  for  which  a  slightly  larger  upper 
bound  is  valid,  by  modifying  some  of  the  good  algorithms. 

3.Sobolev  spaces  and  other  tools.  In  this  Section,  we  will  develop  the  technical  tools 
required  to  prove  the  almost  optimality  of  the  algorithms  to  be  introduced  in  Section  4. 
Many  of  our  tools  can  be  found  in  the  books  by  Grisvard  [17]  and  Necas  [19],  but  there  are 
also  some  results  for  finite  element  spaces  that  do  not  have  continuous  counterparts.  The 
two  books  just  mentioned  consider  Sobolev  spaces  and  related  elliptic  problems  on  non- 
smooth  domains,  normally  Lipschitz  regions.  For  the  problems  in  the  plane  considered  in  this 
paper  the  curvilinear  polygons  considered  in  Grisvard's  book  is  a  natural  family  of  regions. 
These  are  regions  with  corners  with  smooth  curves  connecting  the  corners.  We  will  also 
consider  a  special  class  of  regions  with   boundaries  which  are  not  even  continuous. 

In  most  instances,  we  will  consider  Sobolev  semi-norms  and  norms  defined  on  an 
individual  substructure  or  part  or  all  of  its  boundary.  Without  restricting  the  generality  of  our 
theory,  we  will  principally  consider  the  case  of  second  order  equations  and  Sobolev  spaces 
appropriate  for  that  case.  The  norms  will  contain  certain  scale  factors  which  are  obtained  by 
a  simple  scaling  of  the  standard  norms  defined  on  a  region  with  unit  diameter. 

The  space  ff'  (fl,)  has  a  well  known  subspace  Hq(CI,)  defined  as  the  closure  of  Cq{C1,) 
with  respect  to  the  //'(fl,)  norm.  This  is  a  true  subspace  of  W'(^,)  ^^'^  ^^  'S  also  the  maximal 
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subspace  of  //'(fl,)  for  which  the  extension  by  zero  to  the  complement  of  fl,  defines  a 
bounded  operator  into  /f'(/?^)  .  The  space  of  traces  of  //'(fl,)  is  H"'(dil,)  where  for  an 
arbitrary  curve  F  the  corresponding  semi-norm  and  norm  are  given  by 

K-^n  =  //  l«(^(^))  -  "U(^'))F/k(^')  -  x(s)\^  ds  ds' ;  (3.2) 

r  r 

\\u\\l^,,  =  K^:,  ^  (l/^(r))||u||^^^. 

We  note  that  we  have  a  scale  factor  d(T)  which  is  the  diameter  of  the  curve  T .  The  closure 
of  Cq  (r)  in  the  H'*  (T)  norm  equals  H^(r),  which  turns  out  to  be  equal  to  H^T),  while  the 
maximal  subspace  of  H\r)  for  which  the  extension  by  zero,  to  a  curve  CT  which  is  the 
complement  of  f,  defines  a  bounded  operator  into  H^'^iT  U  CT)  is  known  as  H^q^T)  .  This 
space,  previously  briefly  mentioned  in  Section  2,  is  defined  in  terms  of  a  norm  given  by, 

e 
ll"II^S,(n  =  Kyn  ^  /{l«(^(^))p/^  +  \u{x{s))\yie-s)}ds 

Here  €  =  ((T)  is  the  length  of  the  curve. 

It  is  shown  in  Grisvard  [17],  that  the  last  two  terms  of  this  norm  cannot  be  bounded  in 
terms  of  ||"l|^i/2  p  •    Therefore  Woo(r)  is  not  identical  to  //"^(F).  The  origin  of  the  extra  terms 

in  the  norm  is  explained  in  Grisvard  [17]  and  Necas  [19].  What  is  required  is  a 
straightforward  estimation  of  the  //''(T  U  CT)  norm,  as  defined  in  (3.2),  of  a  function 
defined  on  F  and  extended  by  zero  to  CF.  We  will  carry  out  a  quite  similar  computation  in 
the  proof  of  Lemma  3.2. 

It  is  well  known  that,  the  norm  of  the  quotient  space  h\Q,^)IP (_^,  i.e. 

^m/^JIu-PlU,  (3.4) 

is  equivalent  to  the  //'(fl,)  semi-norm  given  for  €  =  1  in  formula  (3.2).  Here  Pj^  is  the  space 
of  polynomials  of  degree  k.  We  note  that  this  is  a  standard  tool  used  in  error  analysis  for 
finite  elements;  cf.  Ciarlet  [8]. 

We  note,  that  //^(F)  can  also  be  defined  by  using  the  K-method  of  interpolation. 
Thus, 

H^n  =  [//i(F),L2(F)]^, 

cf.  Lions  and  Magenes  [18,  p.  66]. 

The  following  lemma,  plays  an  important  role  in  the  theory. 
Lemma  3.1.   Let  u^  6  V''(fl,)   be   a  finite  element  function  defined  on  a  substructure  of 
diameter  H  and  let  the  diameter  of  the  smallest  element  be  h.  Then  there  exists  a  constant  C, 
independent  of  u^,  such  that 

|2         <  rc\  4-  inatuiuww,.  ii2 


l"'.llnn,)  ^  <^(1  +  log{H/h))\\u,\\l^^^^ 


This  result  is  often  given  for  piecewise  linear  finite  elements  only,  but  it  holds  for  much  more 
general  cases  too;  cf.  Bramble  [5],  Bramble,  Pasciak  and  Schatz  [7],  Oganesyan  and 
Rukhovets  [20],  Thom^e  [22]  or  Yserentant  [26]. 
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As  we  already  know  from  our  discussioD  in  Section  2,  it  is  important  to  be  able  to 
estimate  the  strain  energy  associated  with  a  discrete  harmonic  function  on  a  substructure  in 
terms  of  the  strain  energy  of  its  neighbors.  Consider  a  union  of  four  triangles  where  the 
middle  triangle  fl,  shares  an  edge  with  each  of  its  three  neighbors.  Let  u;,  be  a  continuous 
finite  element  function  which  is  discrete  harmonic  in  H,.  We  note  that  the  union  of  the  three 
outer  triangles  is  not  a  Lipschitz  region;  cf.  Grisvard  [17]  or  Necas  [19]  for  a  definition  of 
Lipschitz  boundaries.  The  following  result,  which  essentially  is  an  extension  theorem  and 
does  not  have  a  continuous  counterpart,  holds. 

Lemma  3.2  The  strain  energy  of  the  triangle,  measured  by  the  square  of  the  //'(fl,)  semi- 
norm,  exceeds  the  sum  of  the  strain  energies  of  the  three  neighboring  triangles  by  at  most  a 
f actor  C  (1 +;og(ff /A  ))^. 

We  first  observe  that 

We  note  that  we  are  considering  only  one  substructure  at  this  time  and  this  result  therefore 
follows  from  the  extension  theorem  discussed  above;  cf.  Bjdrstad  and  Widlund  [5]  or 
Widlund  [24]  for  details.  We  now  use  the  definition  of  the  trace  norm  and  find  that,  since 
an,  =  UV,j, 

Consider  one  of  the  terms  with  ji'k  and  use  some  elementary  inequalities.  This  term  can  be 
estimated  by 

^ti    o 

^11  0  *iy     0 

5  2//    K(x(^))  -  u,{x{Q))M^-s'Ydsds-  +  2/   /    K(Ar(5'))  -u^{x{0))Hs- s'fdsds' 


«-<,. 


'«  0 


s  2  /K(x(5))  -  u,{x{m^/sds  +  2    /  |u,(x(.'))  -  u^{x{QWs'ds'. 

Divide  each  resulting  integral  into  two  and  use  Lemma  3.1.  Here  h,j  is  the  distance  to  the 
next  mesh  point. 

/|«,(X(5))    -    «,(x(0))P/5  d,    =    /|U,(X(5))    -    U,(x(0))P/.d5 
0  0 

+  ^\u,{x{s))  -  u.ixm^/s  ds  ^ 


"u 


C{l+logiH/h)\\u,i)  -  «,(x(0))||^ 


(n,) 
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=  Cil+  log(H/h))\u,i)^  p^-  (u,ix{0))  -  po)  il^.,^^ 
^C(l-logiH/h)r-\Uf,+  PoWl^^^^ 

The  proof  of  the  lemma  is  now  completed  by  a  simple  quotient  space  argument. 

We  will  also  need  to  study  certain  finite  element  functions  that  vanish  at  the  vertices  of 
the  substructures. 

Lemma  3.3.  Let  C  be  linear  operator  in  V''(n,)  which  is  bounded  in  L*(n,).  Assume 
further  that 

Cp  =  0,  p  ^  Po, 
and  that  Cu^  vanishes  at  the  vertices  of  fl,.    Then 

and 

Proof:  The  proof  of  the  first  inequality  follows  exactly  the  same  pattern  as  that  of  Lemma 
3.2  and  no  details  are  therefore  given.  Once  we  recall  the  definition  of  the  Hqq  norm,  we  see 
that  the  second  inequality  requires  very  little  additional  work. 

We  have  mostly  worked  with  semi-norms  in  this  section.  Results  such  as  Lemma  3.2 
can  be  extended  to  the  case  when  the  energy  is  measured  in  terms  of  |«Alivn^  "*"  H^'iHi^nA  > 
by  proving  a  variant  for  the  scaled  norms  introduced  in  (3.2)  and  then  taking  a  convex 
combination  of  the  resulting  inequality  and  that  of  Lemma  3.2. 

4.  Some  almost  optimal  algorithms  for  problems  with  many  substructures.  In  this  section 
we  will  develop  a  number  of  algorithms,  and  prove  the  following  result. 

M&ln  Theorem  The  condition  number  of  every  algorithm  developed  in  this  section  is 
bounded  by  C(l  +log{H/h)^. 

We  begin  by  considering  a  direct  extension  of  the  Neumann-  Dirichlet  algorithm 
introduce  in  Section  2.  We  assume  that  we  can  color  the  substructures  by  two  colors,  red  and 
black.  This  means  that  no  pair  of  elements  which  have  a  common  edge  have  the  same  color. 
Such  a  coloring  is  of  course  not  always  available. 

Denote  the  union  of  the  red  and  black  substructures  by  fl'"  and  Q'-'  respectively.  We 
now  apply  the  Neumann-Dirichlet  algorithm  introduced  in  Section  2  to  these  sets.  When  we 
solve  a  Neumann  problem  on  ft^'\  we  determine  values  of  the  degrees  of  freedom  at  the 
edges  of  all  the  substructures  including  the  vertices  of  the  substructures.  It  is  therefore  easy 
to  see  that  information  will  be  transported  globally  in  this  step  which  involves  a  region  from 
which  roughly  every  second  substructure  has  been  eliminated.  Once  this  Neumann  problem 
has  been  solved,  Dirichlet  problems  are  solved  on  all  the  black  substructures  using  the  traces 
of  the  solution  of  the  Neumann  problem  as  boundary  data.  These  Dirichlet  problems  are 
independent    of    each    other    and    they    can    be    solved    simultaneously    by    using    different 


-11 


processors,  if  available. 

As  in  the  simple  case  of  two  substructures,  discussed  in  Section  2,  the  rate  of 
convergence  is  determined  by, 

The  proof  of  the  almost  optimality  follows  by  applying  Lemma  3.2.  to  estimate  the 
strain  energy  of  each  black  substructure  in  terms  of  its  red  neighbors  and  by  adding  the 
result. 

There  are  several  economical  ways  of  implementing  this  algorithm.  We  first  assume 
that  the  Neumann  problem  on  each  substructure  is  positive  definite.  The  linear  system 
corresponding  to  the  Neumann  problem  on  fl"^  has  about  half  as  many  variables  as  the 
original  problem  on  ft  .  However,  because  of  the  very  special  form  of  ft'",  we  can  decouple 
the  local  problems  in  an  inexpensive  preprocessing  step.  We  can  rewrite  the  Neumann 
problem  on  ft*"  as, 

inf  2      {an ("a."/,)  -  J  f'*h<^^  -  Us  "a  ^^}' 

redn,         '  n,  n, 

subject  to  Uf,  being  continuous  at  all  vertices  of  the  coarse  mesh. 

We  account  for  the  constraints  by  introducing  one  Lagrange  multiplier  for  each 
unknown  at  the  vertices  of  the  substructures.  When  we  write  out  the  resulting  linear  equation 
to  determine  the  values  of  the  nodal  values  of  the  finite  element  function  and  the  Lagrange 
multipliers,  we  obtain  an  indefinite  symmetric  system.  The  coefficient  matrix  has  a  large 
principal  minor  corresponding  to  all  the  nodal  values  and  this  matrix  is  a  direct  sum  of 
matrices  which  correspond  to  finite  element  problems  on  individual  red  substructures.  Each 
column  corresponding  to  a  Lagrange  multiplier  has  only  two  non-zeros.  We  can  therefore 
proceed  to  compute  the  Schur  complement  corresponding  to  the  Lagrange  multipliers.  This 
Schur  complement  is  a  sparse  matrix,  corresponding  to  a  seven  point  stencil  in  a  standard 
case,  and  can  be  factored  into  its  Cholesky  factors  at  an  acceptable  expense  unless  the 
number  of  substructures  is  very  large.  Once  this  work  has  been  carried  out,  an  iteration  step 
involves  the  use  of  these  factors  in  a  back  solving  step.  We  also  need  to  solve  one  problem 
on  each  substructure.  Details,  including  operation  counts,  will  be  given  in  Dryja, 
Proskurowski  and  Widlund  [16]. 

As  an  alternative,  we  can  solve  this  large  system  by  block  Gaussian  elimination.  It  is 
easy  to  see  that  wc  get  a  very  similar  method  if  we  order  the  unknowns  corresponding  to  the 
vertices  of  the  substructures  last.  As  in  the  other  method,  we  have  a  very  large  principal 
minor  which  is  the  direct  sum  of  matrices  corresponding  to  individual  substructures.  When 
using  this  method,  we  need  only  assume  that  a^^^^{uf,,u^)  is  positive  definite. 

Wc  next  consider  a  family  of  methods  which  includes  one  developed  by  Bramble, 
Pasciak  and  Schatz  [7]  and  a  special  multigrid  method  due  to  Yserentant  [26].  Inspired  by 
Yserentant,  we  partition  the  finite  element  functions  into  low  and  high  frequency 
components. 

«A  =  fn^h  +  ("A  -  fH>*h)-  ('*-2) 
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Here  Iff  is  an  interpolation  operator  associated  with  piecewise  linear  continuous  elements  on 
the  coarse  mesh  defined  by  the  substructures.  The  underlying  finite  element  problem  can  be 
defined  in  terms  of  a  more  complicated  finite  element  space. 

Remark.  By  choosing  the  simplest  coarse  mesh  interpolant  the  resulting  algorithm  will  be  as 
simple  as  possible.  By  examining  our  proof  and  the  assumptions  of  Lemma  3.3,  it  will 
become  clear  how  to  extend  the  study  to  a  more  general  choice  of  Iff  . 

We  assume  that  Uf,  is  a  piecewise,  discrete  harmonic  function.  It  is  easy  to  show  that 
IffU^  and  Uj^-IffUf,  belong  to  the  same  class.  Since  the  second  term  in  (4.1)  vanishes  at  all 
vertices  of  the  coarse  mesh,  it  represents  higher  frequencies;  the  pitch  of  a  drum  increases 
substantially  if  held  fixed  at  a  number  of  points.  A  preliminary  preconditioner  is  introduced 
in  terms  of  the  bilinear  form 

flnCV;..  h^'h)  +  <^n  ("a"  '//"a-  "h  -Ih  "O  •  (4-2) 

The  resulting  linear  system  of  equations  is  block  diagonal  if  we  work  with  the  standard  basis 
functions  for  the  piecewise  linear  functions  on  the  coarse  mesh  and  a  complementary  space 
built  by  standard  basis  functions  of  V*".  We  note  that  for  the  basis  functions  of  the  first  kind 
the  second  term  in  (4.2)  vanishes  while  the  first  vanishes  for  those  of  the  second  kind.  The 
first  block  corresponds  to  the  original  problem  solved  on  the  coarse  mesh.  If  this  system  is 
solved  directly,  then  it  provides  the  necessary  quite  modest  amount  of  global  transportation 
of  information.  The  second  term  should  be  changed  since,  by  itself,  it  gives  rise  to  a  linear 
system  which  is  almost  as  difficult  as  the  original  problem. 

As  we  have  noted  the  first  term  in  (4.2)  assures  us  of  global  transportation  of 
information.  If  we  would  replace  the  matrix  representing  the  first  term  in  (4.2)  by  an  identity 
matrix,  destroying  this  feature,  the  condition  number  of  the  resulting  iteration  matrix  would 
increase  by  a  factor  VH^.  Technically,  the  use  of  a  method  where  the  interaction  between 
the  substructures  is  local  only  would  make   it  impossible   to   remove   terms  of  the  form 

Before  we  replace  the  second  term  of  (4.2),  it  is  important  to  note  that  for  Uf,  =  Vf^  the 
expression  in  (4.2)  is  bounded  from  below; 

This  follows  from  the  triangle  inequality.  An  upper  bound  for  the  first  term  is  given  by 

M^I«Hn,  ^C(l^/.g(/f//.))iu|^^„_^  (4.4) 

which  is  obtained  by  using  Lemma  3.1.  and  an  elementary  estimate.  By  using  the  triangle 
inequality,  the  second  term  of  the  right  hand  side  of  (4.3)  is  bounded  similarly.  We  can 
therefore  conclude  that  the  preconditioner  defined  by  (4.2)  has  a  condition  number  on  the 
order  of   (1+  log(H/h)). 

We  will  now  explore  several  different  ways  of  replacing  the  second  term  in  the 
quadratic  form  (4.3).  To  arrive  at  the  algorithm  of  Bramble,  Pasciak  and  Schatz  [7]  this 
second    term    is    replaced    by    terms    related    to    the    trace    of    the    function,    resulting    in   a 
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preconditioner  defined  by, 

'^dhH'h^h)  +  2(«A  -  Va/a  -  VA)«^(r  )  (4-5) 

The  proof  of  the  standard  bound  for  the  condition  number  follows  directly  from  Lemma  3.3. 

In  order  to  derive  the  algorithm  due  to  Yserentant,  we  recursively  split  the  second  term 
in  (4.1), 


where 


^."a  =  'H'*h 


l\»h  =  haH 


and  finally 


Ij^h 


The  preconditioner  is  now  obtained  by  replacing  the  second  term  in  (4.2)  by 

k=\  xi  J/j\  Jl/j.i 

Here  M^  is  the  set  of  mesh  points  after  k  successive  refinements  of  the  original  mesh  defined 
by  the  substructures.  Thus  A/t\.V^_[  is  the  set  of  mesh  points  introduced  in  the  kth  step  of  the 
refinement.  For  a  discussion  of  the  triangulations  for  which  this  argument  applies  and  a  proof 
that  this  algorithm  has  a  condition  number  on  the  order  of  C{\+ log{H/h)y-,  see  Yserentant 
[26]. 

A  hybrid  algorithms  can  now  be  derived.    Combining  the  ideas  of  the  two  previous 
algorithms  we  obtain  a  preconditioner  of  the  following  form, 

e 


+  SK-/f«A.VA-V,)^^^^^^^. 


Here  rf^  correspond  to  a  division  of  Cl  into  finer  substructures.  It  presents  no  further 
difficulties  to  establish  that  the  condition  number  of  this  algorithm  is  on  the  order  of 
(l  +  log(H/h))^.  This  algorithm  might  have  certain  advantages  over  the  two  from  which  it 
originates.  It  might  be  advantageous  to  decrease  the  size  of  the  substructures,  in  particular  if 
the  resulting  linear  systems  of  equations  are  solved  by  a  method  with  a  work  estimate  that 
grows  faster  than  linear  in  the  number  of  unknowns  especially  if  it  can  be  done  without 
increasing  the  cost  of  the  coarse  mesh  computation.  Domain  decomposition  algorithms  might 
also  offer  advantages  compared  with  multi-grid  algorithms  in  that  the  granularity  is  coarser 
and  that  the  algorithm  therefore  is  likely  to  require  less  synchronization. 
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It  is  also  natural  to  replace  the  second  term  in  (4.2)  by  the  quadratic  form  of  the 
Neumann-Diricfalet  algorithm.  Since 

it  is  appropriate  to  introduce  a  scale  factor.  The  resulting  preconditioner  is 

It  is  easy  to  show,  by  using  Lemma  3.2.,  that  the  iteration  operator  has  a  condition  number 
which  is  bounded  by  C(l+  log(H/h))^.  This  algorithm  is  quite  interesting  since  it  does  not 
require  the  use  of  the  J  operator  introduced  in  Section  2  and  used  by  the  algorithm 
introduced  by  Bramble,  Pasciak  and  Schatz  [7].  Compared  with  the  Neumann-Dirichlet 
algorithm  discussed  in  this  Section  it  has  the  advantage  that  the  problems  on  the 
substructures,  represented  by  the  second  term  of  (4.7)  are  uncoupled.  This  follows  from  the 
fact  that  u;,-IffUf,  vanishes  at  the  vertices  of  the  coarse  mesh.  Consideration  such  as  these 
are  even  more  important  in  the  three  dimensional  case,  where  the  vertices  of  a  plane  problem 
correspond  to  edges  with  a  substantial  number  of  points.  We  note  that  the  three  dimensional 
case  differs  considerably  from  the  two  dimensional  problems  considered  here.  Thus  Lemma 
3.1  is  far  from  being  true  in  that  case;  cf.  Dryja  [12]. 

We  conclude  this  paper  by  showing  how  these  algorithms  can  be  extended  to  the 
biharmonic  Dirichlet  problem.  We  first  note  that  the  Dirichlet  conditions  essentially  amount 
to  giving  values  for  the  components  of  the  gradient  Vu^,  on  the  boundary.  The  variational 
formulation  is  given  in  the  Sobolev  space  H^(,il)  and  by  the  trace  theorem  of 
Vuf,  €  (//"'(an))^  .    We  have  the  following  inequality,  for  biharmonic  functions. 

The  extension  of  the  Neumann-Dirichlet  algorithm  and  a  proof  of  almost  optimality  is 
immediate.  When  turning  to  the  other  algorithms,  it  is  relatively  easy  to  understand  what 
properties  are  required  of  the  coarse  mesh  interpolant  IffU^.  We  insist  that  Vlff  is  bounded  in 
L''(n,),  and  that  Iffiax  +  by  +  c)  =  ax  +  by  +  c  for  all  choices  of  a,  b  and  c.  We  can  then 
use  Lemma  3.3.  An  example  of  a  method  which  satisfies  these  requirements  is  given  by  the 
18  degrees  of  freedom  Bell  triangle.  This  conforming  element  is  discussed  in  Ciarlet  [8].  To 
construct  I^  we  use  the  same  finite  elements  on  the  coarse  mesh  defined  by  the  substructures 
and  interpolate  the  values  of  u^  and  Vm;,  while  setting  all  the  second  derivatives  of  Ifju,,  to 
zero.  Certain  bounds  on  the  related  basis  functions  have  to  be  established  in  order  to  prove 
the  L''  estimate.  This  is  done  in  a  computation  resembling  the  study  of  Argyris'  triangle 
given  in  Ciarlet  [8].   Otherwise  the  proof  proceeds  just  as  in  the  second  order  case. 
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